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Abstract. The instability of Bondi-Hoyle-Lyttleton accretion, observed in numerical simulations, is analyzed through known 
physical mechanisms and possible numerical artefacts. The mechanisms of the longitudinal and transverse instabilities, estab- 
lished within the accretion line model, are clarified. They cannot account for the instability of BHL accretion at moderate Mach 
number when the pressure forces within the shock cone are taken into account. The advective-acoustic instability is considered 
in the context of BHL accretion when the shock is detached from the accretor. This mechanism naturally explains the stability 
of the flow when the shock is weak, and the instability when the accretor is small. In particular, it is a robust proof of the 
instability of 3D accretion when y = 5/3 if the accretor is small enough, even for moderate shock strength (M ~ 3). The 
numerical artefacts that may be present in existing numerical simulations are reviewed, with particular attention paid to the 
advection of entropy/vorticity perturbations and the artificial acoustic feedback from the accretor boundary condition. Several 
numerical tests are proposed to test these mechanisms. 
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1. Introduction 

The phenomena described by numerical simulations are 
usually highly simplified compared to the physical reality. 
However simplified, these phenomena are sometimes compli- 
cated enough to challenge our ability to understand them. The 
flow of gas onto a gravitating accretor moving supersonically 
is a classic astrophysical problem (Hoyle & Lyttleton 1939, 
Bondi & Hoyle 1944) which is relevant in many astrophysical 
contexts such as wind fed X-ray binaries, supermassive black 
holes, star formation, and also galaxies in a cluster (see the re- 
cent review by Edgar 2004). Early numerical simulations in 2D 
(Matsuda et al. 1987, Fryxell & Taam 1988) revealed that this 
flow is unstable. 3D simulations (Matsuda et al. 1991, Ruffert 
& Arnett 1994) confirmed the unstable character of the flow, 
displaying however a weaker variability than in 2D. The in- 
stability is not understood even in the simple case of an ideal 
uniform gas. Is this instability physical, or numerical? Is it the 
same instability mechanism in 2D and 3D? How can this unsta- 
ble behaviour be extrapolated to small accretor sizes, which are 
relevant for wind fed X-ray binaries but out of reach of numer- 
ical simulations? Some authors recently doubted that this in- 
stability is physical (Pogorelov, Ohsugi & Matsuda 2000). The 
most recent relativistic simulations by Font & Ibanez (1998a, 
1998b) and Font, Ibanez & Papadopoulos (1999) also showed 
stable flows. What about the published instability mechanisms 



which have been proposed over the years? What is left of the 
longitudinal and transverse instabilities of the accretion line 
(Cowie 1977, Soker 1990, 1991), the "relatively simple" mech- 
anism based on the shock opening angle by Livio et al. (1991), 
the "vortex shedding in the Von Karman manner" (Koide et al. 
1991, Matsuda et al. 1991, 1992)? The present uncertain situa- 
tion reveals how unconvincing the proposed mechanisms were. 
Some were inconclusive, lacked quantitative criteria, and oth- 
ers might also have been incompletely understood. The present 
work aims at clarifying the different instability mechanisms 
and confront them with existing numerical simulations, pay- 
ing particular attention to possible numerical artefacts. 
The general trends based on the existing simulations, and the 
proposed mechanisms are recalled in Sect. 2. The longitudi- 
nal and transverse instabilities of the accretion line are revis- 
ited in Sect. 3 and Sect. 4. The advective-acoustic instability is 
adapted to the BHL flow in Sect. 5. This enables a new look at 
the simulations in Sect. 6, in an attempt to reconcile them. The 
basis of new simulations, free of numerical artefacts, testing 
these ideas, is described in Sect. 7. For the sake of the clarity 
of the paper, the main text contains only the most important 
equations, which summarize the analytical arguments proven 
in appendices A to H. 
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Table 1. Overview of the published numerical simulations of BHL accretion using a polytropic equation of state and a totally ab- 
sorbing accretor: plane accretion in 2D, 3D axisymmetric accretion and full 3D accretion. The first column contains an abreviated 
reference which can be found in the bibliography. 



ref. 


grid 


Mach 


accretor 


index 


transverse 


stability 






number 


r,/r A 


r 


gradients 




3D 


R99 


Cart. 


1.4 - 10 


0.02 - 1 


1.01, 4/3, 5/3 


den 


no 


R97 


Cart. 


0.6 - 10 


0.02 - 1 


4/3,5/3 


vel 


no 


R96 


Cart. 


0.6 - 10 


0.02 - 1 


1.01 




no 


R95 


Cart. 


0.6 - 10 


0.02 - 1 


4/3 




no 


RA95 


Cart. 


3 


0.01 - 10 


5/3 


vel 


no 


RA94 


Cart. 


3 


0.01 - 10 


5/3 




no 


R94 


Cart. 


0.6 - 10 


0.02 - 10 


5/3 




no 


IMS93 


cyl. 


3 


0.125 


isothermal 


den/vel 


no 


MIS92 


Cart. 


3 


0.1 


1.005 -5/3 




no 


MSS91 


Cart. 


3 


0.06 - 0.25 


5/3 




no 


SMA89 


curv. 


1.4 


0.1 


5/3 


vel 


quasi 


LSK86 


Cart. 


3, 16 


0.15 


7/6-5/3 


den 


quasi 


SLK86 


Cart. 


2,4 


0.15 


1 


den 


yes 


3D axisym. 


POM00 


polar 


3- 10 


0.05 


1.01, 1.4,5/3 




yes 


FI98a 


polar 


0.6 - 10 


0.1-2.4 


1.1,4/3,5/3 




yes 


KMS91 


polar 


1.4-10 


0.005 - 0.015 


5/3 




no 


MSS89 


polar 


1.4 


0.01 - 0.05 


5/3 




no 


SMA89 


curv. 


1.4 


0.1 


5/3 




yes 


PSS89 


polar 


0.6-5 


0.125 


1.1,4/3,5/3,2 




yes 


FTM87 


polar 


1.4-4 


0.016-0.13 


5/3 




no 


SMT85 


polar 


0.6-5 


0.1 


1.1,4/3,5/3 




yes 


H79 


polar 


0.6 - 3.6 


0.01 


4/3 






H71 


polar 


0.6 - 2.4 


0.01 


5/3 






2D planar 


POM00 


polar 


3- 10 


0.05 


1.4,5/3 


-/den/vel 


yes 


POM00 


polar 


3- 10 


0.05 


1.01 




no/yes 


POM00 


polar 


4 


0.05 


4/3 


den 


no 


FIP99 


polar 


5 


0.25 


4/3,5/3,2 




yes 


FI98b 


polar 


3-10 


0.25 


1.1,4/3,5/3 




yes 


SMA98 


polar 


1 - 16 


0.005 - 0.05 


isothermal 




no 


BLT97 


polar 


4 


0.001, 0.005 


4/3 


-/den/vel 


no 


ZWN95 


Cart. 


3 


0.03-0.13 


5/3 




yes 


ZWN95 


Cart. 


4 


0.03-0.13 


4/3 


den 


no 


BA94 


SPH 


3 


0.04-0.13 


1.1, 1.3, 1.5 




no 


IMS93 


Cart. 


3 


0.13 


isothermal 


den/vel 


no 


MIS92 


Cart. 


1.4-10 


0.04 - 0.3 


1.005 -5/3 




no 


MSS91 


Cart. 


3-5 


0.06 - 0.25 


1.2,5/3 




no 


SMA89 


curv. 


3 


0.1 


1.5,2 


vel 


no 


TF89 


polar 


4 


0.037 


4/3 


vel 


no 


FT88 


polar 


4 


0.037 


4/3 


den 


no 


MIS 87 


curv. 


1-5 


0.03 - 0.6 


4/3,5/3 


binary 


no 


ABM87 


SPH 


3 


0.13,0.15 


1.5 


vel, den 


yes 



2. An overview of existing simulations of BHL 
accretion and proposed instability mechanisms 

2.1. Numerical simulations are numerous 

Numerical simulations of the BHL problem started with the 
work of Hunt (1971). The instability first appeared in Matsuda 
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et al. (1987), Fryxell & Taam (1988). The many subsequent 
simulations are listed in Tabled Simulations involving more 
complicated ingredients such as realistic heating and cooling 
(e.g. Blondin et al. 1990, Taam, Fu & Fryxell 1991) are not in- 
cluded for the sake of simplicity. The simulations of Table ^ 
are divided into three groups corresponding to plane accretion, 
axisymmetric accretion and full 3D accretion. In each of these 
groups, the listing in chronological order follows the progress 
in computing speed and numerical techniques over the last 
30 years. This progress enabled the simulation of smaller and 
smaller accretors, improving the first attempts by a factor 10. 
Denoting by r A = IGMjv 1 ^ the accretion radius of an accretor 
of mass M and velocity v M , the most recent simulations reach 
rJr A = 0.005 in axisymmetric flows (KMS91), r.Jr A = 0.001 
in planar accretion (BLT97), and r„/r A = 0.01 in 3D accretion 
(RA94). Global trends can be summarized as follows: 

-The shock is always attached to the accretor in simulations 
of 2D planar flow. By contrast, 3D simulations revealed a de- 
tached bow shock, ahead of the accretor, if the accretor is small 
enough and y > 4/3. A calculation in Appendix A suggests 
that the shock should be detached in planar flows with y ~ 3. 

-Although the strength of the instability varies from one 
code to another, the instability is found in numerical simu- 
lations using any of the coordinate systems, Cartesian, polar, 
cylindrical or special curvilinear, even with SPH. If numerical, 
the phenomenon is not specific to a particular grid or a specific 
method. 

-Simulations of plane accretion exhibit the most unstable 
behaviour. The shock moves sideways in a flip-flop manner 
(MIS87). The instability is strongest for small accretors, pos- 
sibly for intermediate Mach numbers (Ai = 3 according to 
MIS87 and POM00). 

-The presence of velocity or density gradients in the trans- 
verse direction of the flow was considered in the earliest unsta- 
ble simulations (MIS87, FT88, SMA89), but MSS91 realized 
that this ingredient is not crucial for instability. This conclusion 
was challenged by ZWN95 who confirmed the flip-flop insta- 
bility when the accretor is a square single cell (as in MSS91), 
but found a stable flow when the accretor is spatially resolved 
and modelized as a polygon. From their point of view, the in- 
stability is related to the detachement of the shock. 

-Stable planar accretion flows were also found by FI98b 
and FIP99 who considered a rather large ratio r„/rA = 0.25 
and POM00 whose method is further discussed below. 

-Axisymmetric flows are generally stable, with few excep- 
tions. Among them, SMT85, FTM87 and MSS89 showed vor- 
tex shedding when the accretor is a hard, non absorbing sphere 
: in this case the problem is similar to the classical flow around 
a sphere, modified by gravity. An important exception to the 
stability of axisymmetric flow is KMS91, who also found vor- 
tex shedding for an absorbing accretor, if Ai > 2.4. It can be 
noted that the accretor size they considered is the smallest ever 
used in axisymmetric simulations. 

-The instability of full 3D accretion is never as strong as 
the flip-flop observed for planar flows, but seems present in 
all published simulations, at least for small enough accretors. 
Even the earliest simulations of LSK86 and SMA89 showed 



some persistant oscillations. The instability seems to be more 
violent if the shock is detached (y = 4/3 and 5/3). 

2.2. Several mechanisms were proposed 

The question of the stability of BHL accretion could in princi- 
ple be solved by performing a perturbation analysis on a sta- 
tionary solution such as obtained by POM00. This procedure 
would be very heavy, and has never been achieved. Using var- 
ious simplifications, six different physical instability mecha- 
nisms have been proposed so far. In chronological order: 

(i) longitudinal instability of the accretion line (Cowie 
1977), 

(ii) transverse instability of the accretion line (Soker 1990, 
1991), 

(iii) shock opening angle (Livio et al. 1991), 

(iv) vortex shedding in the Von Karman manner (Koide et 
al. 1991,Matsudaetal. 1991, 1992), 

(v) local Rayleigh-Taylor (RT) and Kelvin-Helmhotz (KH) 
instabilities (Foglizzo & Ruffert 1999), 

(vi) advective-acoustic cycle (Foglizzo & Tagger 2000, 
Foglizzo 2001, 2002). 

Some of the proposed mechanisms (iii, iv, v) are interesting 
ideas, which need to be developed, but which are not conclu- 
sive at present: 

-The calculation of Livio et al. (1991) is not a stability anal- 
ysis of the shock surface, but rather an attempt to express in 
equations the idea that the pressure should decrease along the 
shock surface. Since no growth rate or typical timescale is com- 
puted, extrapolating on the possibility that this is responsible 
for the violent and chaotic behaviour observed in 2D simula- 
tions is not convincing. 

-Vortex shedding, by the interaction of the incoming gas 
with the "atmosphere" captured by the accretor, was demon- 
strated by SMT85 and MSS89 in simulations where the accre- 
tor is non-absorbing. The fate of this instability in BHL accre- 
tion, where the gas is absorbed at supersonic velocity, is rather 
speculative. 

-The local analysis of FR99 investigates two natural causes 
of instability (RT and KH), and concluded that these are not 
quantitatively convincing without a feedback mechanism. 

Among the cited mechanisms, the only conclusive stabil- 
ity analysis are those of Cowie (1977) and Soker (1990, 1991) 
in the approximation of the accretion line model. This sim- 
plification is known to be valid only when the shock open- 
ing angle is very narrow, i.e. at very high Mach number. Even 
then, it misses the possible interaction between a bow shock 
and the accretor (as described by the advective-acoustic cycle). 
However simplified, these instability mechanisms are physical. 
They should be understood well enough to predict what is left 
of them beyond the accretion line model. 
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Table 2. Accretion line models in 2D and 3D. 




3. A new look at the longitudinal instability of the 
accretion line 

3.1. Physical cause of the longitudinal instability 

The mechanism of the longitudinal instability of the accretion 
line was briefly explained by Cowie (1977) as being due to the 
effect of accreted momentum on density perturbations. Soker 
(1990) challenged this explanation by assessing that this in- 
stability mechanism is independent of accretion and is a mere 
consequence of the acceleration of the flow. A closer look at 
the equations, in Appendix B, shows the weakness of this argu- 
ment. Let us consider more generally the following dynamical 
system: 



dp 

at 



dpv 
dv 



dv 

dt + V 'dr 



H(r), 
F(r, v,p), 



(1) 
(2) 



where H, F are regular functions describing the mass input and 
the force per unit mass. The particular case of the accretion line 
model is described in Tableland Appendix A, where velocities 
are in units of v M , distances r are in units of the accretion radius 
ta , and the line density p is normalized using the mass accre- 
tion rate (the normalization of distances and densities used by 
Cowie (1977) are different by a factor 2). A linearization of 
Eqs. (I1I2> gives the differential equation satisfied by a pertur- 
bation of the mass flux h = podv + voSp (Appendix B). In what 
follows, the subscript for unperturbed quantities Vn,po is omit- 
ted. The solution is written at high frequency u> using the WKB 
approximation : 



P 

v^ 

V V dp ) 



exp 



rt dF dF\ dr 
J \lh ~ P ~dp~)2~v i 

fdr \-i i ClpdFV- 

J v±ir ' jWp) dr 



(3) 



The flow is thus unstable at high frequency if the force per unit 
mass F, acting on the accretion line, depends on density. The 
instability does not depend on the accretion of mass, as stressed 
by Soker (1990), in the sense that it does not depend on the 
function H. Nevertheless it does depend on the accreted mo- 
mentum through the function F. In this sense, accretion plays 
a crucial role in this instability, as initially sketched by Cowie 
(1977). Contrary to the conclusions of Soker (1990), accelera- 
tion within the accretion line is not crucial for this instability 
(see a counter example in Appendix B). In the 3D accretion 



Fig. 1. Schematic view of the accretion line model. Stream lines 
are drawn as solid lines with arrows. When pressure is taken 
into account (Sect. 13. 21 . density perturbations in the accretion 
line propagate as acoustic waves (wavy lines with arrows). 



line model, Eq. Q becomes: 



h(r) oc 



log 4 r 
h(r) oc r ! exp 



■ exp 



1 + i i 3 
iur + (2w) 2 log 3 r 



for r » a, 



2 3 l-i(8aj\~- 



for — - « r « (?, 



(4) 



(5) 



where a is the distance of the stagnation point. The amplitude 
of high frequency perturbations increases far from the accretor 
(r » a). Perturbations in the accreting part of the flow (r < 
a) are also unstable at high frequency co down to a point r oc 
6lT? where the amplication ceases and the WKB approximation 
breaks down. In a numerical simulation, the treatment of high 
frequency perturbations is limited by the numerical resolution. 
According to Eqs. \4l5l . the better the resolution, the stronger 
the instability on both sides of the stagnation point. Formulae 
computed in Appendix B for 2D flows show slight differences 
which are not significant on the scale of a few accretion radii. 

3.2. The longitudinal instability modified by pressure 
forces: analogy with radiation driven winds 

The region of the stagnation point is necessarily subsonic, 
and longitudinal pressure forces cannot be neglected there. 
Pressure forces in the stationary accretion line were considered 
by Wolfson (1977), Yabushita (1978a) and Horedt (2000). The 
simplest formulation corresponds to the isothermal hypothesis, 
in which the dynamical equations dH2l are changed into: 



at dr 

dv dv c 2 dp 

— + v — = F(r,p,v) — , 

dt dr p dr 



(6) 
(7) 
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where c is the isothermal sound speed. Again, a single differ- 
ential equation is obtained for the radial structure of the lin- 
ear perturbation h of the mass flux (Eq. ( IC.41 ). describing the 
propagation of acoustic waves modified by the external forces 
F. The solution is approximated at high frequency through a 
WKB analysis, away from the sonic points (v = +c): 



h ~ p 2 exp 



C dr T. 1 (dF pdF\ 
J v ±c 2\dv c dp J 



for a> » 



(8) 



pdF) 
c dp, 
pv dF 

c 1 dp 

This stability analysis resembles that of radiation driven winds, 
studied by Mestel, Moore & Perry (1976) and Mathews (1976). 
Considering a uniform gravitational acceleration g and a lin- 
earized radiative force (F = Ap - g, H = 0), they showed 
that acoustic waves propagating outwards are amplified by ra- 
diation. In addition to the density effect dF/dp, the velocity 
effect dF/dv is formally well known in line-driven winds (see 
e.g. Carlberg 1980, and more recent reviews by Owocki 1994, 
Feldmeier & Owocki 1998). dF/dv < in the accretion line 
model, implying that this velocity effect is always stabilizing, 
independently of the direction of propagation. This could have 
been anticipated directly from the Euler equation, since a pos- 
itive perturbation of velocity results in a decreased external 
force. By contrast, the effect of the density dependence (term 
dF/dp) is opposite for outgoing and ingoing waves: 
dF pdF 
~dv~ + 



dp 



't l± v ^Zl\ in 3D, 
a \ c I 



(9) 



in 2D. 



(10) 



ri (r2 — a) 

The density effect dF/dp < is thus stabilizing for waves prop- 
agating outwards and destabilizing for waves propagating in- 
wards. This can be understood as follows: with dF/dp < 0, 
a positive density perturbation is associated with a decreased 
external force. According to the Euler equation, this decreased 
force has a damping effect on the positive velocity perturba- 
tion associated with a wave propagating outwards, whereas it 
amplifies the negative velocity perturbation associated with a 
wave propagating inwards. In contrast with the instability with- 
out pressure found by Cowie (1977), this possible amplifica- 
tion of ingoing acoustic waves is not oscillatory. Altogether, 
the density and velocity effects damp outgoing acoustic waves. 
According to Eqs. © and (191 101 . acoustic waves propagating 
inwards may be amplified only if 

c + v<v O0 . (11) 

This condition cannot be fulfilled far from the accretor since 
v ~ Voo for r » a. From this we conclude that the only possible 
amplification of high frequency acoustic waves is restricted to 
ingoing waves in a region of finite size. The size of this region is 
independent of frequency. The amplification factor J\, deduced 
from the WKB analysis, is also independent of the perturbation 
as long as its frequency is high enough to satisfy Eq. ©. It can 
be estimated as follows in 3D: 



J{ ~ exp 



dr 



2c 



for a) » 



r — a v — c 



r - a 



< e 



At 



(12) 
(13) 



This contrasts with the instability found by Cowie (1977) 
which could be arbitrarily fast at high frequency (Eas. l4l5l >. 

3.3. The longitudinal instability beyond the accretion 
line model 

The longitudinal instability of the accretion line proves that if 
the Mach number is high enough, the amplification of ingo- 
ing acoustic waves should be visible in numerical simulations. 
Such a transient amplification, however, cannot be considered 
to be a convincing mechanism to explain the instability ob- 
served in numerical simulations for moderate Mach numbers 
Aloo ~ 3 -5 . A true instability would require a feedback loop, in 
order to build an acoustic cycle. Ingoing acoustic waves may be 
partially reflected outwards near the accretor. Outgoing waves, 
however, are likely to escape to infinity rather than be reflected 
inwards again. Moreover, their amplitude is damped by the ef- 
fect of the accreted momentum. In conclusion, the longitudinal 
instability of the accretion line does not explain the instability 
of BHL accretion. 



4. A new look at the transverse instability of the 
accretion line 

4.1. High frequency approximation of the transverse 
instability 

Soker (1990) extended the stability analysis of Cowie (1977) 
to the case of transverse perturbations in 2D planar flows. The 
position of the accretion line is described in polar coordinates 
9 = 0(r). The angle ¥ between the tangent to the accretion 
line and the symmetry axis, and the transverse velocity vg are 
related to as follows: 



30 

*F = + r — , 

dr 

ve = r [j t +v d-r r 



(14) 
(15) 



As remarked by Soker (1990), the transverse instability is de- 
coupled from the longitudinal one. The angle of the accre- 
tion line satisfies a second order differential equation, which is 
approximated at high frequency in Appendix D using a WKB 
analysis: 



1 r 1 1 n 

°c — exp \iojr + 2 2 (l + i)u> 2 r 2 \ for 

r'i 

jdr 2(1-0 



© oc — exp 

r» 



1 5 



5a 2 



r » a, 



for r <K a. 



(16) 
(17) 



According to Eq. ilbt . the larger the distance from the accretor, 
the larger the amplification of perturbations. As for the longi- 
tudinal instability, the amplification in the region of accretion 
ceases close to the accretor at r oc ar 2 l 5 , and thus depends 
on the numerical resolution. These high frequency estimates 
could be used to compare the radial profiles of the transverse 
and longitudinal instabilities of the accretion line in the linear 
regime in 2D flows: the differences between Eqs. dl6ll7l i and 
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Eqs. (IB.7IB.8> are not significant on the scale of few accre- 
tion radii, and cannot be responsible for the dominant longitu- 
dinal instability observed in the non linear simulations of Soker 
(1991). 



4.2. Comparison of the transverse instability with the 
instability of a flag 

More generally, a system satisfying a transverse equation of the 
form 



(18) 



with arbitrary coefficients A, B, C, is unstable at high frequency. 
The differential equation ( ID.3I satisfied by is written in 
Appendix D. If v + 0, the WKB approximation at high fre- 
quency is as follows: 



0OC (i) 4eXP J 



ioj C B lio)B\i 

— + ^r + t~2 ± ( — 3~) 

v 2v Lv L \ v J / 



dr, 



for to » 



Av Slog B 
rB^ dr v 2 



(19) 



This formulation outlines the role of the restoring force B^ in 
driving the high frequency instability. B = -l/(rip) in the 2D 
accretion line model. This mechanism is reminiscent of the in- 
stability of a flag as described by Argentina et al. (2004), where 
the hydrodynamical force acting on the flag is also proportional 
to the inclination *F. The equation describing the transverse mo- 
tion of a flag with infinite flexibility corresponds to the same 
Eq. (Q3J, with v = 0, A = 0, B = -aU Q and C = a, where U 
is the wind velocity and the coefficient a > characterizes the 
aerodynamic force acting on the flag. The solution of Eq. Jl 8i 
when v = is unstable at high frequency if B < 0: 



oc - exp 
r 



C 

1<JL> — 

B 



B 



A 

7b 



dr. 



(20) 



The presence of finite flexibility ("flexural rigidity") in a re- 
alistic flag material would set an upper bound to unstable fre- 
quencies. The comparison with the instability of a flag suggests 
the possible existence of a transverse instability of the subsonic 
region of the stagnation point (v ~ 0). 

4.3. The transverse instability beyond the accretion 
line model 

4.3.1 . Transient growth of the transverse instability in a 
shock cone 

The accretion line model relies on the asumption that the half 
angle 6q of the shock cone is small compared to other distances. 
Identifying the shock cone with the Mach cone, 6q ~ 1/Aix> is 
indeed small if the incident flow is highly supersonic. With a 
longitudinal wavelength comparable to ~ 2nv / u>, the transverse 
instability of the accretion line model is valid only for 



2nv 2r 
» 



(21) 



Fig. 2. Transverse instability growing in the vicinity of the stag- 
nation point. The shock (thick line), the sonic surface (dashed 
line) and the stream lines (thin lines with arrows) are drawn. 
The thin dotted line delineates the accreted gas. The transverse 
motion of the shock on one side is transmitted to the other side 
through acoustic waves (wavy lines with arrows). 



The maximum exponential amplification 31 deduced from 
Eqs. ( I16I17> is thus bounded by: 



m < exp(2nM ca )'- . 



(22) 



Taking into account the finite width of the shock cone sets an 
upper bound oc Aloo to the frequency of the most unstable per- 
turbations in a plane 2D flow. An additional limitation to the 
accretion line model comes from the acoustic time across the 
shock cone, which is a lower bound to the growth time of the 
instability. This also favours the subsonic region surrounding 
the stagnation point rather than the supersonic regions away 
from it. The instability as described by Soker (1990) is a local 
mechanism, taking place in the advected flow. A global mode 
as observed in the 2D plane simulations of the flip-flop insta- 
bility requires an acoustic feedback inside the subsonic region. 
This instability could be described as a purely acoustic cycle 
between the opposite sides of the shock cone as in Fig.|3] 

4.3.2. Towards a global transverse instability in BHL 
accretion ? 

The amplification of transverse perturbations should be visible 
if the Mach number is high enough. Whether this is sufficient to 
explain the instability observed in 2D flows at low Mach num- 
ber is not clear since this instability is transient. When pressure 
is taken into account, the region of the stagnation point seems 
to be a privileged place for the growth of a global mode involv- 
ing transverse displacement and acoustic propagation. 
The fate of the transverse instability in 3D BHL accretion 
is rather uncertain, as remarked by Soker (1991). Obviously 
transverse motions of the accretion line are forbidden in 3D, 
since incoming symmetrical trajectories which do not intersect 
the displaced accretion line meet along the symmetry axis and 
generate a new accretion line. This does not exclude a possi- 
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ble unstable oscillation of the shock cone in a transverse di- 
rection, such that the symmetry axis stays inside the accretion 
shock. According to Eqs. dl8H9> . the existence of a transverse 
restoring force proportional to the inclination angle is enough 
to generate a high frequency instability. Beyond geometrical 
factors, this force is still present in the 3D geometry. This insta- 
bility mechanism could thus be present for high enough Mach 
numbers, although it has never been observed clearly in 3D nu- 
merical simulations. A correct description of the stability with 
respect to transverse oscillations would require taking into ac- 
count pressure effects and the 3D dynamical deformation of the 
shock cone, which is much more difficult than the accretion line 
formalism. 



5. The advective-acoustic mechanism in BHL 
accretion 

5.1. Schematic formulation of a global cycle 

The advective-acoustic instability deals with the cycle of ad- 
vected perturbations (entropy, vorticity) coupled to acoustic 
waves, between the shock and the sonic surfaces. The coupling 
at the shock is a local process associated to the conservation 
laws through the shock. By contrast, the coupling due to the 
inhomogeneity of the subsonic flow occurs all the way from 
the shock to the sonic surface. F01 showed in a simple radial 
geometry that this acoustic feedback is described by an inte- 
gral over the subsonic region, dominated by the region close to 
the sonic point, where the temperature is highest. This allows 
us to decompose the advective-acoustic cycle in four steps as 
follows: 

(1) advection of an entropy /vorticity perturbation from the 
shock to the sonic point, 

(2) excitation of an acoustic feedback due to the inhomo- 
geneity of the flow, 

(3) propagation of this acoustic feedback towards the 
shock, 

(4) excitation of a new entropy/vorticity perturbation on the 
shock surface. 

Each of these steps j = 1 to 4 is characterized by an efficiency 
Qj measuring the amplification of perturbations. This decom- 
position is motivated by the existence of invariants, allowing a 
direct calculation of Qi based on the conservation of entropy, or 
Qj based on the conservation of acoustic energy. The stability 
of the global cycle depends on the product Q: 



Q = Q\ x Q 2 x <3 3 x <2 4 



(23) 



The linear growth rate of the advective-acoustic cycle, mea- 
sured by the imaginary part of the eigenfrequency u> = (u> r , w,), 
can then be approximated by 



cot ~ — log |<3|, 

T 



(24) 



where t is the duration of the advective-acoustic cycle, gener- 
ally dominated by the advection time. Note that this schema- 
tized approach neglects the purely acoustic cycle, which can 
influence the stability threshold (FTOO, F02). 



5.2. Efficiencies <3, based on radial accretion 

In a radial flow with y > 1, we choose to measure, in the 
entropic-acoustic cycle, the amplification of the perturbation 
/ of the Bernoulli constant : 



/ = vSv + 



7-1 



-cdc. 



(25) 



Let us denote by f ± its values for pressure perturbations 8p ± 
corresponding to an energy flux F ± propagating with (+) or 
against (-) the stream, and f s for entropy perturbations 6S . At 
high frequency, for low degree waves / = 0, 1, 2: 



F oc 



1 



rl/1 2 , 



Ale 2 1 

/* ~ (1+M)c 2 ^, 

yp 

c 2 

f s ~ -SS. 

y 



(26) 
(27) 
(28) 



Non radial entropy perturbations 6S are associated to vortic- 
ity perturbations Sw in shocked spherical accretion through 
Eqs. JE. 17l > to ( IE. 19L so that the amplification of the entropy 
5S — > 5p~ — > 6S' also measures the simultaneous amplifi- 
cation of vorticity 5w — > 5p~ — » 5w' . In the isothermal limit 
(y — » 1), vorticity is more appropriate than entropy to de- 
scribe the advective-acoustic cycle, which becomes a vortical- 
acoustic cycle (F02) . 

The steps (1) and (3) of advection and propagation in the 
spherical accretion flow studied by F01, F02 are deduced from 
Eqs. \26\ and d28l > and the conservations of entropy SS and 
acoustic energy F~: 



fS 2 

J son L son f L son 

for Wr „ 

Jsh C sh rson 

f- 7Vl sh ' 
J son L son 



(29) 
(30) 



where the subscripts "son" and "sh" refer to the sonic point and 
the shock respectively. Although the acoustic energy is con- 
served, the amplitude of outgoing waves f decreases ((33 < 1) 
due to the geometric dilution in a diverging flow. 
The advection of an entropy perturbation in a hot region may 
greatly increase the thermal energy carried by this perturbation, 
by a factor comparable to the temperature ratio (Eq. (I29l >). As 
sketched by FTOO, the difference of energy is carried away by 
acoustic waves, propagating both upward and downward. Even 
the waves propagating downward may be refracted upward at 
low enough frequency (F01). The ratio ./son//son is mus of order 
unity: 



s \ _ /on i 

Q 2 = -f ~ 1. 

J son 



(31) 



The advective-acoustic coupling Q4 at the shock is deduced 
from a local analysis of a perturbed shock, performed in 
Appendix F: 



f s 

£sh 
/sh 



1 - At 



sh 



At 



(32) 
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This equation is identical to Eqs. (16) and (C.ll) in FTOO, ob- 
tained for radial perturbations. Q4 may reach large values in 
nearly isothermal flows with strong shocks (A1 s h ~ 1/VMi), in 
which case the dominant advected perturbations are non radial 
(F02): this is the basis of the vortical-acoustic cycle. The fac- 
tor (1 - Alsh) in Eq. J32t is responsible for the damping of the 
advective-acoustic coupling for weak shocks. 
Altogether, 



Q • 



Cson 1 - M 



sh 



Csh 



Ml 



(33) 



This approximate formula illustrates the two regimes of effi- 
cient advective-acoustic coupling identified by FTOO, F01, F02: 

(i) Strong temperature gradients are responsible for an effi- 
cient triggering of acoustic waves from advected entropy per- 
turbations, which is the basis of the entropic-acoustic cycle 
(FTOO, F01). 

(ii) Strong entropy/vorticity perturbations can be produced 
at the shock if the post shock Mach number A4 s h is small. 
These two sources of amplification can be used as guidelines 
for anticipating the properties of advective-acoustic cycles in 
BHL accretion flows. However, the dependence of Q on the fre- 
quency and the degree / of the perturbation requires further cal- 
culations. In particular, F01 showed that the acoustic feedback 
is strongest for non radial modes I = 1 . Note also that Eq. J33l > 
is singular if c son — > 00 as is the case for 3D accretion with 
7 = 5/3, which deserves a more careful analysis (Appendix G). 

5.3. From radial to BHL accretion 

5.3.1 . Extrapolation to detached shocks 

The locus of the advective-acoustic cycle is different if the 
shock is attached or detached (Fig.[3Jl. A global mode between 
a bow shock and the accretor resembles the advective-acoustic 
cycle in a shocked radial flow. In this case the efficiency of 
the advective-acoustic cycle can be estimated from the study of 
radial accretion. By contrast, if the shock is attached to the ac- 
cretor, most of the flow is accreted from behind in a supersonic 
manner. The value of Q2 in this particular geometry cannot be 
directlly extrapolated from its value in a radial flow. 

5.3.2. Geometrical factors 

In BHL accretion, the shape of the shock surface is not only 
non-spherical but open to infinity. The value of Q3 might be 
reduced by a geometrical factor ~ 2 compared to spherical ac- 
cretion because a significant fraction of acoustic waves may 
propagate away from the region of accretion and leave the cy- 
cle (Fig.|5Jl. 

Conversely, the value of Qi should be increased by the ampli- 
fication of vorticity perturbations through the local KH and RT 
mechanisms (FR99), due to vorticity and entropy gradients in 
the post-shock flow. 

The efficiencies <3, estimated in spherical geometry should thus 
be considered, at best, as very rough approximations of those 
in the BHL flow. 





Fig. 3. Schematic view of a global mode if the shock is de- 
tached or attached. Drawn lines have the same meaning as in 
Fig- El Perturbations of entropy/vorticity (circular arrows) fol- 
low the flow lines from the shock to the accretor, producing an 
acoustic feedback which propagates to the shock. A fraction of 
the acoustic energy may leak outside of the accretion cylinder. 



5.4. Influence of the parameters 7, vH» and r* 

According to FTOO, F01, the stability of the entropic-acoustic 
cycle depends essentially on the temperature increase between 
the shock and the sonic point. The most unstable cycle involves 
high frequency acoustic waves, those able to explore the hottest 
parts of the flow but still be refracted out, with a wavelength 
slightly larger than the smallest size of the sonic surface. Using 
the Bernoulli equation, the temperature on a point of the sonic 



surface is directly related to its distance r son to the accretor: 



y-i 

y+l 



2GM 



+ ct \M 



7- 



1 



(34) 



The closer the sonic surface to the accretor, the higher the tem- 
perature, and the more efficient the entropic-acoustic cycle. The 
calculation of Appendix A indicates that the sonic surface is al- 
ways attached to the accretor if 7 = y max , with 



7n 



= 3 in 2D, 



7max = g in 3D, 



(35) 
(36) 
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even if the shock is detached. Fig. 17 of R94 shows the shape 
of the attached sonic surface for 3D accretion with y = 5/3. 
Together with Eq. j33> . the strength Q of the instability should 
increase when y approaches y max and r» decreases. The strength 
of the instability should be asymptotically independent of the 
incident Mach number for strong shocks. Conversely, the insta- 
bility should be suppressed if the shock is weak (A1 s h ~ 1). 
According to Appendix A, a critical index y cr ; t < y max must 
exist, below which the shock is always attached to the accretor, 
whatever its size. Numerical simulations with y — 4/3 (R95) 
suggest that y cl ; t < 4/3 in 3D flows. The entropic-acoustic cy- 
cle is expected to be an efficient instability mechanism in the 
range [y c rit.7max]> as l° n g as the distance of the sonic surface is 
small enough. 

Nearly isothermal flows (y ~ 1) could be unstable through 
the vortical-acoustic cycle, fed by Q4 » 1 for strong shocks. 
However, the effect of the acoustic feedback in the particu- 
lar geometry of an attached shock is uncertain. Neither the 
vortical-acoustic mechanism nor the extrapolated transverse 
instability manage to explain why isothermal BHL accretion 
seems so much more unstable in planar flows (SMA98) than in 
3D simulations (R96). 



6. Can the different numerical simulations be 
reconciled? 

6.1. Some simulations are stable 

The stability of the 2D planar flow simulated by FI98b, FI99, 
with a relativistic accretor contrasts with the many unsta- 
ble simulations performed in Newtonian gravity. This is by 
no mean a relativistic stabilization of BHL accretion through 
relativistic effect, as recognized by FI98b. Indeed, they re- 
stricted their studies to Vc/ciight = 0.5, which corresponds to 
a rather big Schwarzschild radius in units of the accretion ra- 
dius: rsch/'A = 0.25. These flows would have been stable even 
in the Newtonian limit. 

According to ZWN95, the flow simulated by MSS91 becomes 
stable if the square accretor is replaced by a polygon. Although 
the shape of the accretor may influence the instability thresh- 
old, the existence of strongly unstable simulations in polar co- 
ordinates (BLT97 and SMA98) shows that 2D accretion can be 
unstable even if the accretor is perfectly spherical. The smaller 
accretor size and higher resolution used by BLT97 and SMA98, 
compared to ZWN95 (see Table [0, may be a hint in favour of 
an advective-acoustic mechanism. A direct comparison, how- 
ever, is hampered by the fact that the adiabatic indices are dif- 
ferent in these three simulations. The influence of the shape 
of the accretor, demonstrated by ZWN95, speaks against the 
transverse acoustic instability. These rather indirect arguments, 
if not conclusive, show at least that it is possible to reconcile 
existing simulations of 2D plane accretion in the framework of 
the advective-acoustic mechanism. 

The stable 3D simulations can also be analyzed in the 
framework of the entropic-acoustic cycle. This cycle is stabi- 
lized by a weak shock, and destabilized by a small accretor 
size. 



-Table [2 indicates that most axisymmetric simulations of 
an absorbing accretor are stable (SMT85, PSS89, SMA89), 
with the exception of KMS91 which shows an instability when 
r„/rA ^ 0.05 and M > 2.4. The existence of a threshold for 
the accretor size and the minimum Mach number fits perfectly 
with the entropic-acoustic mechanism. 

-In the axisymmetric simulations of FI98a with Voo/cught 
0.5, the Schwarzschild radius is too big to allow for a detached 
shock. According to Newtonian 3D simulations, the shock dis- 
tance scales like a fraction (e.g. typically 0.2-0.4 in R94) of the 
accretion radius. A slower black hole would have a detached 
shock, and the entropic-acoustic instability could develop nat- 
urally if the temperature gradient is sufficient. The shock gets 
detached in the simulation of FI98a for y = 5/3, v^lc = 0.15, 
but the shock is then too weak (Moo = 1 -5) to be destabilized. 
Indeed, the Newtonian simulations of Ruffert with Moo = 1 .4 
were also stable. A decisive test could be made by simulating 
an axisymmetric flow with y = 5/3, Moo = 3, and v M /c < 0.1. 

-The apparent stability observed in the simulations of 
POM00, in particular for y = 5/3, cannot be explained by sim- 
ple considerations about the accretor size and shock strength. 
This result may be attributed to the particular numerical method 
of local time step used by the authors. In Sect. 3 of POM00, 
the authors make the following statement: "It is important to 
note from the very begining that we seek steady state solution 
and generally do not perform time-accurate calculations (local 
time step inside each cell for the sake of computational effi- 
ciency)". This method might not be adequate to propagate high 
frequency acoustic waves across the subsonic region. 

6.2. Numerical artefacts in the simulations 

Numerical issues are numerous. Besides the damping effect 
of numerical viscosity (SPH and Eulerian codes were com- 
pared by BA94) and the possible axis effect in axisymmetric 
simulations (FTM87), more subtle effects can be understood 
through the advective-acoustic cycle. This instability mecha- 
nism is physical, but may be artificially triggered or damped 
by numerical effects such as the carbuncle phenomenon at the 
shock, the boundary condition at the surface of the accretor and 
the grid size in between. The accuracy of both the advection of 
entropy/vorticity perturbations and the propagation of acoustic 
waves, between the shock and the sonic surface, is crucial for 
advective-acoustic instabilities. 

6.2.1 . Carbuncle phenomenon at the shock 

POM00 drew attention to possible numerical instabilities at the 
shock in the BHL flow. In the region where the shock is par- 
allel to the grid, the carbuncle instability (e.g. Robinet et al. 
2000) may favour the generation of vorticity and entropy per- 
turbations, which in turn can feed an advective-acoustic cycle. 
Conversely, one should carefully check the effect of any numer- 
ical procedure designed to damp the carbuncle instability at the 
shock, since it may also damp the coupling between advected 
and acoustic perturbations. 
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Fig. 4. Numerical damping of a sinusoidal perturbation of en- 
tropy (full line) or vorticity (dotted line) advected over one 
wavelength, in a direction parallel to the grid, as a function 
of the number of grid points per wavelength. The flow is ho- 
mogeneous with y = 5/3, and moves uniformly at Mach 0.2. 
The dashed line corresponds to the damping of an acoustic 
wave propagated over one wavelength. In the three cases, the 
wavevector of the perturbation is parallel to the flow. 



6.2.3. Boundary condition on the accretor 

The size of the accretor is known to play an important role in 
the strength of the instability. One reason for this is that the 
sonic surface ahead of the accretor is generally attached to the 
accretor in numerical simulations. Boundary conditions at the 
surface of the accretor are therefore even more crucial, since 
they are in contact with a subsonic flow. The advection of en- 
tropy and vorticity perturbations through the boundary condi- 
tion should be considered carefully in order to avoid an artifi- 
cial acoustic feedback from the accretor. For example, impos- 
ing a transverse velocity equal to zero in the perturbed Bondi 
flow would generate a spurious feedback f~ from any non 
spherical entropy/vorticity perturbations f s , and artificially ex- 
cite a vortical-acoustic cycle. A calculation in Appendix H of 
this condition at the boundary of a uniform, parallel flow shows 
that 

(4) = -1 if / * 0. (37) 

V J I boundary 

The artificial coupling due to inadequate boundary conditions 
can be as strong as the physical coupling Qi expected from the 
temperature gradients in the flow. 



6.2.2. Grid resolution between the shock and the 
accretor 

Vorticity is not usually computed with as much accuracy as mo- 
mentum or energy in classical numerical schemes such as those 
used for BHL simulations. The artificial generation of vorticity 
at the interface of nested grids could feed a vortical-acoustic 
cycle. 

Conversely, insufficient numerical resolution is responsible for 
an artificial damping of vorticity and entropy waves. As an ex- 
ample, Fig. |4] shows a measure of the damping of the ampli- 
tude of entropy and vorticity perturbations advected over one 
wavelength, in a direction parallel to the grid, as a function of 
the number of grid cells per wavelength. This test is performed 
on a Cartesian grid in 2D using the same PPM technique as 
in Ruffert (1994). The correct advection of entropy and vortic- 
ity in this numerical simulation requires as much as 10-15 grid 
cells per wavelength. Perturbations with a wavelength shorter 
than 10 grid cells are significantly damped over one wave- 
length. The same test performed on acoustic waves propagating 
over one wavelength shows a smaller damping: 5-10 grid cells 
are enough. This is a strong constraint for the correct calcu- 
lation of an advective-acoustic cycle at high frequency which 
involves perturbations with a wavelength comparable to the ac- 
cretor size, when the sonic surface is attached to it. The grid 
size should thus be at least 15-20 times smaller than the ac- 
cretor size in order to properly describe the advective-acoustic 
coupling in the inner regions of the flow. This example illus- 
trates the fact that the advective-acoustic instability can be im- 
peded in a numerical simulation with insufficient resolution. 
Determining how these numbers depend on the numerical tech- 
nique is beyond the scope of the present paper. 



6.3. Physical or numerical instability? 

Many of the numerical artefacts discussed in Sect. l6.2l were not 
considered by the authors of the existing simulations. Could 
new simulations of BHL accretion, corrected from numerical 
artefacts, be stable? The physical arguments of Sect. [5] prove 
that the flow must be unstable, at least in the case y ~ 5/3, 
Moa > 3, for a small enough accretor: in this case, the effi- 
ciency Q of the entropic-acoustic cycle diverges when the ac- 
cretor size r, — > 0. Even the leak of acoustic energy cannot 
diminish the efficiency Q by more than a finite geometrical fac- 
tor 2 - 3, at most. The argument is weaker for y < 5/3, because 
Q depends strongly on the unknown shape of the sonic surface. 

7. Future numerical tests of the instability 
mechanism 

Future simulations should consider carefully the numerical 
artefacts listed in Sect. 16.21 The boundary condition at the sur- 
face of the accretor should be designed to absorb entropy and 
vorticity perturbations silently: this can be tested in a uniform 
flow. An alternative would be to find a set of parameters such 
that the accretor is fully embedded inside the sonic surface. 
Since y — 5/3 and y ~ 1 are ruled out by FR97 in 3D, an in- 
termediate choice could be y — 4/3 with Mca = 3, and a small 
enough accretor. The numerical issue of the absorbing bound- 
ary condition can also be solved naturally with general relativ- 
ity, for any value of y, since the flow is bound to be supersonic 
on the horizon of the black hole (FI98a,b). 

The stable simulations analyzed in Sect. 16.11 call for new 
simulations which would directly test the efficiency of the 
advective-acoustic mechanism in accretion flows where the 
shock is detached: 
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-all the axisymmetric, Newtonian simulations with y — 5/3 
should be unstable if Moo > 3 and a small enough accretor 
( r */ r A ~ 0.005 seems to be enough according to KMS91). 
In particular, new axisymmetric, relativistic simulations sim- 
iliar to FI98a should become unstable when considering a 
slower accretor with a moderate shock, such as Vco Might = 0.1, 
y = 5/3, Mc<, = 3. The axisymmetric simulations of POM00 
should also be unstable if their numerical technique allows for 
the advective-acoustic cycle. 

-2D planar simulation should be unstable through the 
entropic-acoustic cycle if the shock is detached: comparing a 
simulation with y ~ 3 (see Appendix A) and the classic flip- 
flop obtained for y < 5/3 could help to understand the respec- 
tive influences of the entropic-acoustic cycle and the purely 
acoustic transverse instability. 

In an unstable simulation, the advective-acoustic mecha- 
nism could be tested directly by measuring the incoming en- 
tropy/vorticity perturbations and the outgoing acoustic flux in 
the subsonic region of the flow, as measured by Blondin et al. 
(2004) in his simulations of spherical accretion. An indirect 
way of testing the instability mechanism is to measure the ef- 
fect of the various physical parameters. The linear growth rate 
of the entropic-acoustic cycle should increase when the accre- 
tor size is decreased, decrease for a weak shock and saturate for 
strong shocks. The present understanding of the 3D advective- 
acoustic instability when y ~ 5/3 implies that its strength 
for an accretor size r„ ~ 10~ 5 ta may be significantly under- 
estimated by the existing simulations (r„/rA > 10~ 2 ). Tracing 
the power spectrum of the mass accretion rate (e.g. R95), as 
a function of the numerical resolution in the range accessible 
to computers, could give a hint on the extrapolation to smaller 
accretor sizes. 

In order to better understand the instability of isothermal 
flows, it would be interesting to try to discriminate between a 
transverse instability (purely acoustic) and a vortical-acoustic 
cycle. The accretor size and boundary condition should play a 
more important role in a vortical-acoustic cycle than in a trans- 
verse instability. 

8. Conclusion 

The physical mechanisms underlying the longitudinal and 
transverse instabilities of the accretion line have been clarified. 
The destabilizing factor for the longitudinal instability is the 
density dependence of the external force per unit mass. The 
transverse instability is ruled by the restoring force propor- 
tional to the local inclination of the accretion line. The analogy 
with the instability of a flag has been outlined. 
The overstable amplification of longitudinal high frequency 
perturbations of the accretion line found by Cowie (1977) is 
greatly affected by pressure forces. Density perturbations are 
propagated as acoustic waves. Those propagating outwards are 
damped, whereas those propagating inwards are transiently 
amplified. The analogy with the instability of radiation driven 
winds has been drawn. 

The transverse instability of the accretion line is limited if the 
finite width of the shock cone is taken into account. A feedback 
process is necessary to explain the global flip-flop instability 



observed in planar 2D simulations. It seems possible that the 
mechanism of the transverse instability, modified by the propa- 
gation of acoustic waves within the subsonic region of the flow, 
leads to a global unstable mode. However, the influence of the 
shape of the accretor, demonstrated by ZWN95, may be a hint 
against this explanation of the instability. 
For the first time, the advective-acoustic instabilitiy has been 
analyzed in the context of BHL accretion. This analysis is 
based on the extrapolation of the properties established in 
spherically symmetric flows (F01, F02). For this reason, the rel- 
evance of the advective-acoustic instabilities in BHL accretion 
is convincingly demonstrated only when the shock is detached 
from the accretor. The difference of geometry precludes an ac- 
curate prediction of the instability threshold in BHL accretion, 
especially since it is very sensitive to the size and shape of the 
sonic surface. Nevertheless, the analysis is predictive enough 
to assess that the advective-acoustic mechanism 

- must be stable in the limit of a weak shock, 

- must be unstable for 3D accretion flows with y ~ 5/3, 
a reasonable shock strength Moo > 3 and a small enough 
accretor size. 

Numerical artefacts have also been discussed, specifically 
the grid size and the boundary conditions. We comment on 
how these artefacts must be circumvented to produce reliable 
numerical simulations of BHL accretion. 
Surprisingly, it seems that all the existing numerical simula- 
tions can be reconciled in the framework of the advective- 
acoustic instability, with no striking contradiction, even when 
the shock is attached to the accretor. Several numerical tests of 
these ideas have been proposed. 

Besides new numerical simulations, future analytic work 
ought to describe in more detail both the transverse instability 
modified by pressure forces, and the efficiency of the advective- 
acoustic cycle when the shock is attached. This could help un- 
derstand why isothermal accretion is so much more unstable in 
2D than in 3D. 
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Appendix A: BHL stationary flow in 2D and 3D 

A.1. Adiabatic index and pressure forces 

The sonic radius in radial Bondi accretion, deduced from the 
Bernoulli equation and the conservation of mass flux, is differ- 
ent in 2D and 3D: 

3 -yGM 

r so „ = — in2D < (A- 1 ) 
5-3yGM . 

r SO n = —r-^r in3D - ( A - 2 ) 

The existence of a sonic radius in radial accretion on a point 
like accretor requires y < 5/3 in 3D, whereas plane supersonic 
accretion is possible up to y < 3. This illustrates the fact that 
the 3D convergence of flow lines produces stronger pressure 
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gradients than in 2D. These pressure gradients act against grav- 
ity. Extending to plane flows the argument used in FR97 for 3D 
flows, the sonic surface of a stationary BHL accretion must in- 
tersect the sphere of radius rg deduced from Eqs. ( IA.1IA.2> . (w - kv) 2 - — 
with 



rewritten using the stationary flow equation Eqs. (2.7a-2.7b) of 
Soker (1990), as follows: 

2 - /fel "-in2D, (B.2) 



ro = 



(A3) 



This suggests that the shock should be detached in planar ac- 
cretion with y close to 3. SMA89 considered plane accretion 
with y — 2 but the shock was still not detached. The ratio ro/r^ 
is a function of y, M x : 



La 

r,\ 
£0 



7 



M 2 



4 i + l 
5 -3y 



2 

Ml 



Ml 



1 + y -^Ml 



in 2D, 



in 3D. 



(A.4) 



(A.5) 



As already noted in FR97, the sonic surface must also be at- 
tached to the accretor if ro/rA > 1, because the sonic sur- 
face cannot extending beyond the distance ~ r^ of the stag- 
nation point. This concerns in particular isothermal flows with 
a strong shock, for which ro/rA °c M^. A critical index y cr jt 
must therefore exist, below which the shock is always attached 
to the accretor, whatever its size. 

A.2. Accretion line 

The mass flux per unit of length, along the accretion line, is 
constant in 3D whereas it varies like r ~~- for a planar flow 
(Soker 1990). By integrating Eq. Q using Table|2] 



pv = 2(r J - a'-) in 2D, 
pv — r — a in 3D, 



(A.6) 
(A.7) 



where a is the distance of the stagnation point. The asymptotic 
velocity within the accretion line, deduced from Eq. (0 for r » 
a, is different in 2D and 3D: 



v ~ 1-4" in 2D > 

ri 



log r 
1 - — in 3D, 



where A is a constant. 



(A.8) 
(A.9) 



Appendix B: Physical cause of the longitudinal 
instability of the accretion line 

The growth rate computed by Cowie (1977) in 3D is the imag- 
inary part of the complex frequency oj: 



(to - kv) = ik 



dv 1 

dr 2r 2 



The growth rate in 2D is given by strictly the same formula, 
although the accretion terms involve geometrical factors l/ri 
in Table [2] The absence of such factors in Eq. ( IB. Il l led Soker 
(1990) to conclude that this instability is independant of ac- 
cretion. This argument is insufficient, since Eq. dH.lt > can be 



ikl-v . 



2 p 



in 3D. 



(B.3) 



The geometrical factor rJ is then clearly apparent. A closer 
look at the equations shows that the longitudinal instability is 
due to the density dependance of the force acting on the ac- 
cretion line. Using the same normalizations of density, velocity 
and distances as in Soker (1990), the time dependent equations 
correspond to Eqs. dll2> . 

A linearization of Eqs. ill 12b gives the differential equation sat- 
isfied by a perturbation of the mass flux h = pdv + vdp: 



, d 2 h dh 



dr 2 dr 



d v 2 



dF dF 



dr p dv dp 



, f d v dF 

-icohip- to) — —— 

\ dr p dv 



0. 



(B.4) 



The solution is written in Eq. Q at high frequency a) using 
the WKB approximation. Contrary to the conclusions of Soker 
(1990), acceleration within the accretion line is not crucial for 
this instability. The simplest flow in which a similar instability 
occurs would be a flow with uniform density po and velocity 
vo, subject to a force depending linearly on density: F = a(p — 
po) and without any mass accretion (H = 0). The stationary 
flow being uniform, the uniform velocity could also be taken to 
be equal to zero owing to a simple change of reference frame. 
The evolution of perturbations can then be calculated precisely. 
The frequency u> and the wavevector k are related through the 
following dispersion relation: 



(w - kv ) = ikpo 



dF_ 



The growth rate corresponds to the imaginary part w,-: 



OJi = 



apo 



2(to r - kvo) 



(B.5) 



(B.6) 



The amplification of perturbation is thus exponential in the di- 
rection of the external force for a positive enhancement of den- 
sity. In the accretion line model, this amplification is weaker 
due to the weak density dependence of the external force. Using 
the asymptotic behaviour v oc —\ /n close to the accretor, and 
Eq. ( IA.9> in 3D flows leads to Eqs. ( 1415 \ . By contrast, the same 
calculation using Eq. JA.8> in 2D flows leads to: 



i r 1 - 1 1 

h(r) oc r 8 exp \ia>r ± 2(1 + i)a>^A 2 r* , 

for r » a, 



(B.l) ^( r ) oc ri ex P 



2 3 1 - i 1 
— iur 1 + -u) 2 r 

3 2a\ 

1 

for — - <K r <K a. 

a>2 



(B.7) 



(B.8) 



The difference of asymptotic behaviors between the 2D and 
3D cases is not significant. 



T. Foglizzo, P. Galletti, and M. Ruffert: Instability mechanism 



13 



Appendix C: Effect of pressure forces on the 
longitudinal instability 

The linearization of the perturbed flow leads to define the per- 
turbation / of the Bernoulliu constant as follows: 

/ = v6v + c 2 — , 
P 



(CI) 



in order to obtain a simple second order differential system sat- 
isfied by (f,g): 



7 o dh 
(c -v 2 )— = ico(pf-vh), 



dF_ 
3v~ 



dF 
vp-z- 
dp 



(C.2) 



(C3) 



A single differential equation is obtained: 



d 2 h dh 



d 



or 1 or \ or \ p 

dF dF) , f d v dF) n 

+v- p— > + uoh\p- iu> - — } = 0. 

dv dp J { dr p dv 



(C4) 



Before looking for approximate solutions to this equation, the 
effect of pressure forces can be easily incorporated in the sim- 
ple toy model used in Appendix B, with uniform density and 
velocity: 



1 dF , 

to = -— +kv± 

2 dv 



dF 1 IdF 



kc +ikp Tp-A\^ 



(C.5) 



The problem is formally identical to that idealized by Mestel, 
Moore & Perry (1976) and Mathews (1976), for radiation 
driven winds, where F = Ap - g. At high frequency: 



. i (dF pdF\ 
w .* (v±c) + _ _ ± __ . 



(C6) 



The stability of high frequency acoustic waves depends on the 
sign of the quantity between parenthesis. In a non uniform flow, 
a similar conclusion can be reached at high frequency through 
a WKB analysis, leading to Eq. (|8j. 

Appendix D: Transverse instability of the 
accretion line 

The differential equation satisfied by © is: 

d 2 ® d® f <91ogr 2 v 2ia) 1 / I \) 
dr 2 dr 1 dr v r \ vp \ v/f 



© \ 2 n . ■ 1 3 
r < rpco + 2ia>pv + icor 2 ■ 

rpv 2 [ 2n- 
The WKB approximation is thus: 



0. 



(D.l) 







1 +i 
exp + — -co 
2-- 



,v)4 r 
— exp 

i-i j 

1 C dr 

2 I i I I 



ICO 

V 



1 



2r?vp 



dr 



3v 

for u> » — . 

2r 



(D.2) 



The second order differential equation in the general case cor- 
responding to Eq. d 1 81) is: 



2 d 2 ® d® Ivdrv \ 

dr 2 dr \r dr J 

I v A + B 7 \ 
+® \-iu> + iuC - tu r I = 0. 



(D.3) 



Appendix E: Proof that 6K = in the shocked 
Bondi flow 

Let us recall the differential system satisfied by the perturba- 
tions f,g defined in F01 (Eqs. (B18-B19)): 



df iuM 2 f 



iojv 



,6S 



^dg icoM g 



ia)f 



dr l-M 2 c 2 (l-M 2 ) cor 2 



iL 2 i5K 
/+ — 



(E.l) 



(E.2) 



where the constant 5K and the function fi(r, to, I) are defined 
by: 



SK = r 2 v ■ (V x 6w) + 1(1 + l)c 2 — , 

r 

2 _ l(J+l). 2 2, 

V = 1 t-t( c - v )■ 



(E.3) 
(E.4) 



Consider a spherical adiabatic shock with incident Mach num- 
ber Mi in the radial direction. Let this shock be perturbed by 
a sound wave with frequency a> propagating against the flow in 
the subsonic region, producing a displacement A£(6, <p) and a 
perturbation Av(6>, <p) of the radial velocity of the shock. Using 
the index "1" before the shock, and "2" after it, the conserva- 
tion of mass flux and energy across the shock can be written as 
follows: 



P\(v\ - Av) = p2(v2 + 6V2 - Av), 



(E.5) 



(vi - Av) 2 



c\ ( Vl + 5v2-Av) 2 (c 2 + 6c 2 ) 2 

n + T = n + i — ' (K6) 

2 2 T _ l 

where quantities are measured at the position r s h + A£. Keeping 
the first order terms, and using the defnition of /, g, together 
with the entropy equation, we obtain: 



/ = (v 2 ~ Vi)Av, 



— - — |Ar + ().V. 

V 2 V] , 



(E.7) 
(E.8) 



A third equation relating 6S to Av, A£ could be deduced us- 
ing the conservation of impulsion, in the spirit of Nakayama 
(1994). A more direct derivation can be obtained by noting that 
entropy is conserved before and after the shock, and that the 
entropy jump across the shock depends only on the local value 
of the incident Mach number M\ in the frame of the shock 



(Eq. Ell): 



Av 



M[(r sh + AO = Mi(r sh ) + — + A^ 



dMi 



M 



c\ dr 

, s (, irjv 2 \Av 
(r sb )+ 1 + — — , 
V u>r I c\ 



1 = 



Slog M 
Slogr 



(E.9) 
(E.10) 
(E.ll) 
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Using the Rankine-Hugoniot jump conditions, we obtain 



6S = - 



\y(M\ - l) 2 c 2 



\ cor I v\ 



(r + i) 2 M 2 c\ 

Equations (IE.7IE.8> thus become: 



2c 2 . Av 
-i-(M - 1)-, 

' + 1 Vi 



c 2 M 2 



7+lcj 



M, 



1 - 



WV2 



y + 1 cor 



Av 
vi ' 



(E.12) 

(E.13) 
(E.14) 



The perturbation of non radial velocity 5vg, 6v<p immediately 
after the shock is deduced from the conservation of the tangen- 
tial component of the velocity across the shock, in the spirit of 
Landau & Lifschitz (1987). 

Vi - v 2 dA( 



5v, 



r 86 ' 
v\ - v 2 dA{ 



Sv <f> = • „ -i 
r sin 6 dip 



(E.15) 
(E.16) 



The vorticity immediately after the shock is deduced from the 
non radial component of the linearized Euler equations (B. 1 1) 
and (B. 12) of F01, together with Eq. dE~7l and Eqs dE.15IE.16> : 



0. 



d 5S 



Wfl 



rv sin 6 dip y 



rv 80 y 



(E17) 
(E.18) 

(E.19) 



Immediately after the shock, the constant SK defined by 
Eq. JE.3> is computed using Eqs. ( IE. 181 and ( IE. 191 : 



6K = 0. 



(E.20) 



Since 6K is conserved through the Bondi flow, 6K is uniformly 
equal to zero. As a consequence, using the integrated expres- 
sion of the vorticity (Eqs. (B5) to (B7) in F01), the vorticity 
perturbation is described by Eqs. ( IE. 171 to ( IE. 19i throughout 
the flow. 

Appendix F: Advective-acoustic coupling at the 
shock 

The perturbations /, g after the shock are decomposed as fol- 
lows: 



/ = r+r+f 



(F.l) 
(F.2) 



where f s , g s correspond to the entropy/vorticity wave associ- 
ated to the entropy perturbation SS with 6K = 0, and f ± ,g ± 
correspond to the purely acoustic waves propagating in the di- 
rection of the flow (index +) or against the flow (index -). 
Neglecting the coupling between the entropic and acoustic 
waves in the vicinity of the shock, the entropy wave f s , g s is 
advected at the velocity of the fluid: 

df_ 
dr 

(F.4) 
or v 



ICO 



(F.3) 



l CO 



Replacing these derivatives in Eqs. iE. 1IE.2> . we obtain: 
„ 1-M 2 7 6S 



1 -^M\ 



(F.6) 



Acoustic waves are described by Eqs. JE. 1IE.2> in the absence 
of entropy perturbations, i.e. when 5S = 0. Using the WKB ap- 
proximation of F01, the radial derivative of is approximated 
by: 

9/* ICO M + yU 

dr c 1 - M 2 " 
g ± is deduced from Eqs. ( IE. 1> and JF.7I >: 



(F.l) 



Mc 



(F8) 



The linear system iF. 1IF.2I can be transformed using Eqs.( IF.5I 
IF.6I > and !F.8t . in order to express the acoustic perturbations f ±m . 

M 



r 



f± 



-c 2 g-(l±ijM)f s 



(F9) 



Using this equation immediately after the shock, with 
Eqs. dEHl to EH : 



r- 



.y+l c 2 SS M 2 (V 2 + 2Min + Mr 2 ^ 
4y 1 - Mr 2 n 



1 + jlM.2 



From Eqs. JF.5i and JF. 101 . 

f_ = _4 Af il-M 2 ){l-M- 2 ) 

f- y + 1 M 2 (1 - fiM 2 )(fi 2 + 2fiM 2 + M' 2 ) ' 
At high frequency such that fi ~ 1 , 
f_ _ 4 1 +M 2 



1 -Mr 2 



y + 1 1 + 2M 2 + Mr 2 

I -M 2 
M 2 ' 



M 2 



(F10) 
(Ell) 

(F.12) 
(F.13) 



Appendix G: Entropic-acoustic coupling in 3D for 

r~5/3 

The value of |<3| deduced from F01 for y close to 5/3 increases 
with frequency like co 1 / 3 (Eqs. 28-29 of F01), up to a maxi- 
mum reached near the cut-off frequency. This behaviour can 
be understood in the framework of FT00, which argued that 
the efficiency of the entropic-acoustic coupling is related to the 
increase of enthalpy between the shock and the sonic point. 
A slight correction, however, should be made. Rather than the 
naive guess |<3| 2 °c c 2 on /c 2 h of Eq. (23) in FT00, one must take 
into account the fact that the coupling of entropy perturbations 
to acoustic waves must occur before the sonic radius in order to 
allow significant outgoing acoustic flux. The effective radius of 
coupling r e ff was computed analytically in Appendix E of F01 
(Eqs. E6 and El 3): 



\Q\ 2 



c 2 (r eft ) 
c 2 (r sh ) 



(G.l) 
(G.2) 
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The efficiency Q 2 deduced from the analytical calculations of 
F01 indeed scales like the ratio of enthalpies between the shock 
and the effective point of coupling r & g. The radius r e g coincides 
with the wavelength v/a» cut of an entropy perturbation with 
a frequency u> r close to <w cut : the enthalpy effectively "seen" 
by this perturbation is not the enthalpy at the sonic radius but 
rather at the effective radius r er y. From this point of view, using 
the same notations as in F01, \Qs |/=i » \Q$ |/ = o is a natural con- 
sequence of <jj Ut » a>Q Ut : non radial perturbations effectively 
"see" regions of higher enthalpy. The most efficient entropic- 
acoustic coupling is reached at frequencies close to the refrac- 
tion cut-off, where fi ~ 1 . 

Appendix H: Artificial acoustic feedback from the 
boundary condition on the accretor 

The perturbed velocity in the direction perpendicular to the 
flow is deduced from the linearized Euler equations, using the 
vorticity given by Eqs. JE.18IE. 19l : 



A 1 8 f 
6VB = ' ~^Z' 

iu>r oti 

c 1 df 

Sv v = : : — Z7T' 

iu>r sin 9 dip 



(H.l) 
(H.2) 



A combination of these equations involves the eigenvalues of 
the Laplacian in spherical coordinates: 



sinf? 



(H.3) 



Imposing Sv ± — on the boundary consequently requires / = 
there, if the perturbation is not spherically symmetric (I + 0). 
The acoustic feedback f~ associated to the entropy/vorticity 
perturbation f s passing through the boundary condition 6v ± = 
0, is deduced from the decomposition iF.H , with / + = and 
/ = 0: 



= -1 if 1*0. 



(H.4) 



boundary 
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